<?php

    /**
     * @package JAMA
     *    For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
     *    orthogonal matrix Q and an n-by-n upper triangular matrix R so that
     *    A = Q*R.
     *    The QR decompostion always exists, even if the matrix does not have
     *    full rank, so the constructor will never fail.  The primary use of the
     *    QR decomposition is in the least squares solution of nonsquare systems
     *    of simultaneous linear equations.  This will fail if isFullRank()
     *    returns false.
     * @author  Paul Meagher
     * @license PHP v3.0
     * @version 1.1
     */
    class PHPExcel_Shared_JAMA_QRDecomposition {
        const MATRIX_RANK_EXCEPTION = "Can only perform operation on full-rank matrix.";
        /**
         *    Array for internal storage of decomposition.
         * @var array
         */
        private $QR = [];
        /**
         *    Row dimension.
         * @var integer
         */
        private $m;
        /**
         *    Column dimension.
         * @var integer
         */
        private $n;
        /**
         *    Array for internal storage of diagonal of R.
         * @var  array
         */
        private $Rdiag = [];

        /**
         *    QR Decomposition computed by Householder reflections.
         * @param matrix $A Rectangular matrix
         * @return Structure to access R and the Householder vectors and compute Q.
         */
        public function __construct($A) {
            if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
                // Initialize.
                $this->QR = $A->getArrayCopy();
                $this->m  = $A->getRowDimension();
                $this->n  = $A->getColumnDimension();
                // Main loop.
                for ($k = 0; $k < $this->n; ++$k) {
                    // Compute 2-norm of k-th column without under/overflow.
                    $nrm = 0.0;
                    for ($i = $k; $i < $this->m; ++$i) {
                        $nrm = hypo($nrm, $this->QR[$i][$k]);
                    }
                    if ($nrm != 0.0) {
                        // Form k-th Householder vector.
                        if ($this->QR[$k][$k] < 0) {
                            $nrm = -$nrm;
                        }
                        for ($i = $k; $i < $this->m; ++$i) {
                            $this->QR[$i][$k] /= $nrm;
                        }
                        $this->QR[$k][$k] += 1.0;
                        // Apply transformation to remaining columns.
                        for ($j = $k + 1; $j < $this->n; ++$j) {
                            $s = 0.0;
                            for ($i = $k; $i < $this->m; ++$i) {
                                $s += $this->QR[$i][$k] * $this->QR[$i][$j];
                            }
                            $s = -$s / $this->QR[$k][$k];
                            for ($i = $k; $i < $this->m; ++$i) {
                                $this->QR[$i][$j] += $s * $this->QR[$i][$k];
                            }
                        }
                    }
                    $this->Rdiag[$k] = -$nrm;
                }
            } else {
                throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ARGUMENT_TYPE_EXCEPTION);
            }
        }    //    function __construct()

        /**
         *    Return the Householder vectors
         * @return Matrix Lower trapezoidal matrix whose columns define the reflections
         */
        public function getH() {
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    if ($i >= $j) {
                        $H[$i][$j] = $this->QR[$i][$j];
                    } else {
                        $H[$i][$j] = 0.0;
                    }
                }
            }
            return new PHPExcel_Shared_JAMA_Matrix($H);
        }    //    function isFullRank()

        /**
         *    Return the upper triangular factor
         * @return Matrix upper triangular factor
         */
        public function getR() {
            for ($i = 0; $i < $this->n; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    if ($i < $j) {
                        $R[$i][$j] = $this->QR[$i][$j];
                    } elseif ($i == $j) {
                        $R[$i][$j] = $this->Rdiag[$i];
                    } else {
                        $R[$i][$j] = 0.0;
                    }
                }
            }
            return new PHPExcel_Shared_JAMA_Matrix($R);
        }    //    function getH()

        /**
         *    Generate and return the (economy-sized) orthogonal factor
         * @return Matrix orthogonal factor
         */
        public function getQ() {
            for ($k = $this->n - 1; $k >= 0; --$k) {
                for ($i = 0; $i < $this->m; ++$i) {
                    $Q[$i][$k] = 0.0;
                }
                $Q[$k][$k] = 1.0;
                for ($j = $k; $j < $this->n; ++$j) {
                    if ($this->QR[$k][$k] != 0) {
                        $s = 0.0;
                        for ($i = $k; $i < $this->m; ++$i) {
                            $s += $this->QR[$i][$k] * $Q[$i][$j];
                        }
                        $s = -$s / $this->QR[$k][$k];
                        for ($i = $k; $i < $this->m; ++$i) {
                            $Q[$i][$j] += $s * $this->QR[$i][$k];
                        }
                    }
                }
            }
            /*
            for($i = 0; $i < count($Q); ++$i) {
                for($j = 0; $j < count($Q); ++$j) {
                    if (! isset($Q[$i][$j]) ) {
                        $Q[$i][$j] = 0;
                    }
                }
            }
            */
            return new PHPExcel_Shared_JAMA_Matrix($Q);
        }    //    function getR()

        /**
         *    Least squares solution of A*X = B
         * @param Matrix $B A Matrix with as many rows as A and any number of columns.
         * @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
         */
        public function solve($B) {
            if ($B->getRowDimension() == $this->m) {
                if ($this->isFullRank()) {
                    // Copy right hand side
                    $nx = $B->getColumnDimension();
                    $X  = $B->getArrayCopy();
                    // Compute Y = transpose(Q)*B
                    for ($k = 0; $k < $this->n; ++$k) {
                        for ($j = 0; $j < $nx; ++$j) {
                            $s = 0.0;
                            for ($i = $k; $i < $this->m; ++$i) {
                                $s += $this->QR[$i][$k] * $X[$i][$j];
                            }
                            $s = -$s / $this->QR[$k][$k];
                            for ($i = $k; $i < $this->m; ++$i) {
                                $X[$i][$j] += $s * $this->QR[$i][$k];
                            }
                        }
                    }
                    // Solve R*X = Y;
                    for ($k = $this->n - 1; $k >= 0; --$k) {
                        for ($j = 0; $j < $nx; ++$j) {
                            $X[$k][$j] /= $this->Rdiag[$k];
                        }
                        for ($i = 0; $i < $k; ++$i) {
                            for ($j = 0; $j < $nx; ++$j) {
                                $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
                            }
                        }
                    }
                    $X = new PHPExcel_Shared_JAMA_Matrix($X);
                    return ($X->getMatrix(0, $this->n - 1, 0, $nx));
                } else {
                    throw new PHPExcel_Calculation_Exception(self::MATRIX_RANK_EXCEPTION);
                }
            } else {
                throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MATRIX_DIMENSION_EXCEPTION);
            }
        }    //    function getQ()

        /**
         *    Is the matrix full rank?
         * @return boolean true if R, and hence A, has full rank, else false.
         */
        public function isFullRank() {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($this->Rdiag[$j] == 0) {
                    return false;
                }
            }
            return true;
        }
    }
